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^ ; Abstract 

The Rational Hybrid Monte Carlo (RHMC) algorithm extends the Hybrid Monte Carlo algorithm 

> . 

for lattice QCD simulations to situations involving fractional powers of the determinant of the 
quadratic Dirac operator. This avoids the updating increment (dt) dependence of observables 

oo ; 

which plagues the Hybrid Molecular-dynamics (HMD) method. The RHMC algorithm uses rational 



approximations to fractional powers of the quadratic Dirac operator. Such approximations are 
only available when positive upper and lower bounds to the operator's spectrum are known. We 
apply the RHMC algorithm to simulations of 2 theories for which a positive lower spectral bound 
is unknown: lattice QCD with staggered quarks at finite isospin chemical potential and lattice 



QCD with massless staggered quarks and chiral 4-fermion interactions (xQCD). A choice of lower 
bound is made in each case, and the properties of the RHMC simulations these define are studied. 
Justification of our choices of lower bounds is made by comparing measurements with those from 
HMD simulations, and by comparing different choices of lower bounds. 
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I. INTRODUCTION 



For lattice field theories with fermions, and in particular for lattice QCD, where there 
exists a local action describing the dynamics with the desired number of fermion flavours, 
the preferred simulation method is the hybrid Monte-Carlo (HMC) algorithm I]. The HMC 
augments a hybrid molecular-dynamics update with a discrete 'time' step (dt), with a global 
accept /reject step at the end of each trajectory. This accept/reject step makes the algorithm 
exact, removing all dt dependence from measured observables. 

For staggered and improved staggered fermions with numbers of flavours which are not 
multiples of 4 and for other fermion schemes with odd numbers of flavours, there is no 
local action, and one needs to take fractional powers of the fermion determinant. Until re- 
cently, simulations of such theories used hybrid molecular-dynamics simulations with 'noisy' 
fermions, usually the R algorithm (HMD(R)) [2]. The problem is that in HMD(R) simula- 
tions, the observables depend on dt, the leading departure from the desired dt = results 
being 0(dt 2 ). Hence one should either simulate at several (small) dt values and extrapo- 
late to dt = 0, or perform simulations at such a small dt value that the 0(dt 2 ) errors are 

n t;;: R HMc ^ mm — - — _ of ^ 

minant of the quadratic Dirac operator, one takes fractional powers of the operator itself. 
This defines an action (nonlocal) and thus permits a global accept/reject step. Since the 
quadratic Dirac operator is positive definite, such fractional powers are well defined in terms 
of the eigenmodes of the operator. To make this algorithm practical, one needs approxima- 
tions to the fractional powers of the eigenvalues, which can be applied directly to matrices. 
The RHMC algorithm uses rational approximations and their partial-fraction expansions 
to approximate such fractional powers over the range of eigenvalues of the quadratic Dirac 
operator. Provided that the condition number (ratio of maximum to minimum eigenvalues) 
of this operator is not too large, sufficiently good approximations to this fractional power 
can be obtained with a modest order of numerator and denominator polynomials in this 
approximation. 

The first theory we consider is lattice QCD at finite isospin chemical potential ///, but with 
no explicit symmetry breaking parameter ^, |^, Q] . The spectrum of the quadratic Dirac 
operator for this theory has an unknown lower bound. However, we can make a reasonable 
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guess as to what is a conservative estimate of the lower bound for the small fij values of 
interest. The need for such an exact algorithm is indicated, since, in HMD(R) simulations, 
the Binder cumulants which are used to determine the nature of the finite temperature 

nnnn 

phase transition for this theory depend strongly on dt |9|, HO, [11|, 112J . We have adapted the 
RHMC algorithm to this theory with 3-fermion flavours and tested the reversibility needed 
for its implementation. To check the validity of our speculative lower bound to the Dirac 
spectrum, we have compared our measurements with extrapolations to zero dt of our results 
from earlier HMD(R) simulations at various dts. In addition, we have performed simulations 
using rational approximations which assume much smaller lower bounds to the spectrum, 
at the largest /i/ value and selected j3 and m values which we used with the original choice 
of bounds. Indications are that our original choice of lower bounds are adequate. 

The second theory we are considering is 2-flavour lattice QCD at zero quark mass, using 
the xQCD action, which incorporates an irrelevant chiral 4-fermion interaction that allows 
simulations at zero quark mass ^| . Here again, we have reasonable estimates of the up- 
per bound of the spectrum of the quadratic Dirac operator, but the lower bound is unknown. 
Again we make an educated guess as to what is a reasonably conservative lower bound for 
the spectrum of this operator and test our choice. For a 24 3 x 8 lattice, we found that using 
32-bit floating point precision, was the main limit on reversibility. With 64-bit precision, 
the convergence criteria for the multimass Dirac operator inversion was the limiting factor, 
leading us to prefer using 64-bit precision. We tried 3 different choices of speculative lower 
bounds on the quadratic Dirac operator comparing the inversions at the beginning and at 
the end of each trajectory. For the 2 lowest choices, these inverses agreed to within the 
known accuracy of the rational approximations used, for all trajectories in the run. The 
highest choice agreed with the other 2 for almost all trajectories in the run. In addition we 
made direct comparison of the results of simulations using the RHMC algorithm with earlier 
measurements using the HMD(R) algorithm on small (8 3 x 4) lattices. 

Section 2 introduces the 2 theories under consideration with rational approximations. In 
section 3 we discuss the case of QCD at finite /!/, while in section 4 we consider xQCD. 
Finally in section 5 we present discussions and conclusions. 
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II. QCD AT FINITE m, xQCD AND RATIONAL APPROXIMATIONS 



In this section we present two Lattice QCD actions for which the lower bound on the 
spectrum of the quadratic Dirac operator is unknown. This is important for the application 
of the RHMC algorithm for numbers of fermion flavours other than multiples of 8. One 
can only find a rational approximation to a fractional power a of a positive (semi-) definite 
matrix to any desired relative accuracy, if one knows positive lower and upper bounds to its 
spectrum. 

If M is a positive (semi-) definite matrix with eigenvectors |A) and corresponding (non- 
negative) eigenvalues A, and a is a real number, 

M a = £A«|A)(A| (1) 

A 

gives a unique definition of any real power of M. If one knows the spectral bounds on M 
i.e. that < a < A < b, for some a and b, then one can find a rational approximation to X a 
of the form 

A" « P(A)/Q(A) (2) 

where P and Q are polynomials in A, such that for a any given (small) positive e 

Max | A °-P(A)/Q(A)| 
A e K b \ < e - (3) 

Ase— > or a/6 - > 0, the orders of P and/or Q needed — > oo. The method of calculating 
the best approximation for fixed orders of P and Q was specified by Remez. We use an im- 
plementation provided by Clark and Kennedy jl5| (Note that this approximation also yields 
the best rational approximation to A~ a ). This rational approximation to the eigenvalues of 
M a provides a rational approximation to M a itself as 

M° w P(M)Q~ l (M), (4) 

which can be implemented directly in terms of the matrix M without having to find its 
eigenvalues and eigenvectors. Expanding the right-hand-side of equation H] in terms of par- 
tial fractions allows the action of M a on any vector to be calculated using a multishift 



conjugate gradient algorithm 
gradient inversion. 
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171 ] for little more than the cost of a single conjugate 
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A. QCD at finite /z/ 

Lattice QCD at finite isospin chemical potential /zj, has the staggered quark action 

S f=Yl x\Jp(k^m) + m ]x + i^X^2X , (5) 

sites 

where Tpi^r^i) is the standard staggered quark transcription of pwiih the links in the +t 
direction multiplied by exp(|r 3 /z/) and those in the —t direction multiplied by exp(— |t3/Zj). 
For finite temperature simulations with /z/ < m r , we set the coefficient A of the explicit 
symmetry breaking term to zero. As it stands, this action describes 4 quark flavours (in the 
continuum limit). 

To tune this to Nf flavours following the RHMC approach, we replace this with the 
pseudo-fermion action 

S pf = plM- N f/% (6) 
where is the momentum conjugate to the pseudo-fermion field ip. 

M = [lp{\rm) + m] f [^~T3//i) + m] + A 2 (7) 

is the quadratic Dirac operator, and we set A = 0. To apply the RHMC algorithm we will 
need rational approximations to Ai~ Nf ^ 8 and Ai ±N f/ 16 . 

It is easy to see that [cosh (|/zj) + 3 + m] 2 is an upper bound to the spectrum of M.. We, 
in fact, choose an upper bound of 25 which exceeds this bound for any choice of /zj and 
m we are ever likely to consider. At /z/ = 0, m 2 is a lower bound on the spectrum. One 
might expect that this bound is overly conservative, and that the actual lower bound would 
be closer to \m lx . For fij > 0, we do not know the lower bound of the spectrum of the 
quadratic Dirac operator. However, the effective pion mass is m v — fii, so if this gives at 
least some guide as to the lower bound of this operator, we would expect the lower bound of 
its spectrum to decrease smoothly as fij is increased, becoming zero at fij = m^, the value 
of iii at the zero temperature transition to the superfluid phase. Our observations from 
HMD(R) simulations, indicate that in this low /z/ regime the Dirac operator becomes more 
singular as /ij is increased, in a controlled way. Since we are interested in quark masses 
m > 0.02, we choose a speculative lower bound of 1.0 x 10~ 4 , 4 times lower than the known 
lower bound at /zj = for m = 0.02. Tests are applied to see that this is reasonable. We 
restrict ourselves to the case Nf = 3. 



B. xQCD 



The xQCD staggered lattice action incorporates an irrelevant chiral 4-fermion interac- 
tion of the Gross-Neveu form, which allows us to run at zero quark mass. In terms of 
pseudo-fermion fields ijj and their conjugate momenta the fermion action, modified for 
implementation of the RHMC algorithm is 

S P f = P U^A)' N f/% (8) 

where 

A= ^ + m +^E( (J i + ief i) ( 9 ) 

with i running over the 16 sites on the dual lattice neighbouring the site on the normal 
lattice, e = ^—i) x +y+ z + f anc l p j s the usual gauge- covariant "d-slash" for the staggered 
quarks, o and it are the auxiliary fields which implement the 4-fermion interactions. The 
action for these auxiliary fields is 

S^E-A^ + tt 2 ), (10) 

s 8 

where the sum is over the sites of the dual lattice. For more details of the rest of the 
action we refer the reader to our earlier work Q]. We did, however make one minor change, 
multiplying the (a, n) 'kinetic energy' by 7, which acts as a large mass and slows down the 
dynamics of these chiral fields. 

Although we do not know the upper bound of the spectrum of the Dirac operator, be- 
cause of the interaction with the chiral field, we know from our HMD(R) simulations that 
(|A" / 7(<j 2 + it 2 )) is just above 1. For N f = 2 and 7 = 10 the coefficient \Np = 2.5. 
This plus the fact that the a and ir interacting with the fermion fields at a given point are 
averages over 16 points on the dual lattice, suggest that 1 is a conservative upper bound 
for the contribution of the chiral fields to the eigenvalues of the Dirac operator. We thus 
assume an upper bound of 25 for the spectrum of the quadratic Dirac operator. If this were 
violated, it would only affect modes close to the ultraviolet cutoff on the theory, which are 
of limited interest. We shall later present evidence that the magnitude of the chiral field 
which interacts with the fermion field at a point is bounded by 1. 

More important is an accurate knowledge of the lower bound on the spectrum of the 
quadratic Dirac operator. Here we simply assume that the action of the chiral fields on the 
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spectrum is similar to giving the quarks a small mass. Our speculative lower bound is that 
for a regular staggered action with a quark mass 0.001, i.e. 1 x 10~ 6 , and we perform tests 
of this assumption. 



III. TESTS OF THE RHMC FOR LATTICE QCD AT FINITE m. 

The RHMC algorithm for QCD at a finite chemical potential /// for isospin is implemented 



following Clark and Kennedy 
Michael Clark's Remez code 



The rational approximations are obtained using the 
We have been using this method for simulating Nf = 
3 flavour QCD. Here we use a (20,20) rational approximation to .M*- 3 / 16 ** needed at the 
beginning of each trajectory, and a (20,20) rational approximation to A4 <_3//16 ' ) needed at 
the end of the trajectory to calculate the final 'energy' (value of the classical molecular- 
dynamics Hamiltonian describing the system), over the interval [1 x 10~ 4 ,25] discussed in 
section 2. These approximations have a maximum relative error of 6.1 x 10~ 12 . The evolution 
of the fields over the trajectories were performed using a (10, 10) rational approximation to 
_yVj(-3/8) over i n t er val [1 x 10~ 4 ,25], which has maximum relative error of 4.4 x 10~ 6 . 
Since we have been using relatively small lattices 8 3 x 4 and 12 3 x 4, we have found 32-bit 
floating-point precision (accumulations are performed in 64-bit precision) to be adequate. 

For the RHMC algorithm to be exact, i.e. for it to produce an ensemble of configurations 
which have the correct Boltzmann distribution for the given action with no dt dependence, 
the updating algorithm must be reversible. The updating method specified by Kennedy et 
al. is reversible for infinite precision arithmetic, even when the multimass conjugate gradient 
solver used to invert the Dirac operator is stopped prior to convergence. (The inversions at 
the beginning and end of each trajectory still need to be exact for this to hold). With the 
finite precision used by computers, this reversibility must be tested. As mentioned above, 
we use single precision (32-bit) floating point arithmetic for simulations of lattice QCD at 
finite [ij. To test this we performed a 1000 trajectory run on a 12 3 x 4 lattice at (3 = 5.1285, 
m = 0.03 and fii = 0.3, in which every second trajectory was run forward and then reversed, 
returning to its initial state, dt = 0.05 for this experiment, and the trajectory length At = 1. 
\ii = 0.3 is chosen, because it is the largest fii used in our production runs, m = 0.03, just 
above the critical mass, is the smallest mass for which we run at /i/ = 0.3, and f3 = 5.1285, 
close to the transition /3 for these values of m and ///, is one we use for production runs. 
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We expect therefore that the Dirac operator will be as close to singular as we will encounter 
in any of our simulations, and that if the algorithm is satisfactory here, it will work for the 
whole range of parameters we use. 

The change in 'energy' 5E, for each of the 500 reversed trajectories is plotted in figure^,. 
The change in energy for each of the non- reversed trajectories is shown in figure ^p. We 
choose E as a measure of reversibility because it is the quantity that is used in the global 
Metropolis accept reject step. Hence departures of 5E from zero, its value for true reversibil- 
ity, indicate the size of errors introduced from departures from exact reversibility coming 
from finite precision errors potentially magnified by instabilities. As we see, these are indeed 
small in magnitude - < 3 x 10~ 3 , and 3 orders of magnitude smaller than the change AE 
in energy over a typical trajectory used for normal updating. In addition, the average SE 
over the 500 trajectories considered is only « 1.8 x 10~ 4 , indicating that there is no sign of 
a systematic bias, so that the small errors introduced by this lack of reversibility are likely 
to cancel. 

Most of the simulations of lattice QCD at finite isospin density, which we have run using 
this code use dt = 0.05, and all use trajectory length At = 1. With this choice we find an 
acceptance rate of ~ 70% for the generated trajectories. Most of our runs used the same 
(3 for the updating as for the global Metropolis step, since we did not find large gains from 
choosing different /3s. 

We check our speculative lower bound on the spectrum of the quadratic Dirac operator 
of 1 x 10~ 4 , at our largest fii (fii = 0.3) at both masses (0.03 and 0.035) where we performed 
production simulations. This we did by running using rational approximations, valid over 
intervals with the same upper bound, but different lower bounds, and compared observables. 
At m — 0.035, fii = 0.3 and j3 = 5.1370, we performed a 300,000 trajectory run with rational 
approximations valid over the interval [1 x 10 -4 , 25] used for production running and a 
300,000 trajectory run with the same rational approximation for updating, but a (20, 20) 
rational approximation valid over the range [1 x 10 -5 , 25] for initiating each trajectory and 
for calculating the energy at the end of each trajectory. This approximation has maximum 
relative error of 2.0 x 10~ 10 . A comparison of the observables from these 2 simulations 
is given m table UK. At m = 0.03, /i 7 = 0.3, f3 = 5.1290 we ran one 300,000 trajectory 
run with rational approximations valid over our default interval ([1 x 10~ 4 ,25]) and one 
300,000 trajectory run with rational approximations valid over the interval [1 x 10 -6 , 25] 
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12 d x4 lattice, 0=5.1285, m=0.03, ^=0.3 
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12 d x4 lattice, £=5.1285, m=0.03, ^=0.3 
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FIG. 1: a) The energy change for 'closed' trajectories (t — > t + 1 — > i). b) The energy change for 
'open' (t—*t + l) trajectories. 



for initialization, measurement and updating each trajectory. For this later interval we used 



(25, 25) rational approximations with maximum relative error 2.0 x 10~ n at the ends of each 
trajectory and (15, 15) rational approximations with a maximum relative error of 7.2 x 10~ 7 
for the updating. Observables from these 2 runs are compared in table |TJd. 



m = 0.035 


, iii = 0.3, (3 


= 5.1370 


Spectrum range 


[1 x 10" 4 ,25] [1 x 10 -5 ,25] 


B 4 


1.860(36) 


1.879(41) 


Pc 


5.13744(24) 


5.13787(22) 


plaquette 


0.51086(66) 


0.51134(60) 


Wilson Line 


0.3300(74) 


0.3245(68) 




0.6274(94) 


0.6347(86) 


,•3 

Jo 


0.0450(10) 


0.0442(10) 


m = 0.03, 


fjt! = 0.3,p-- 


= 5.1290 


Spectrum range 


[1 x 10~ 4 ,25] [1 x 10~ 6 ,25] 




1.704(30) 


1.722(39) 


Pc 


5.12826(18) 


5.12820(21) 


plaquette 


0.50776(60) 


0.50753(68) 


Wilson Line 


0.3773(68) 


0.3792(76) 


(H>) 


0.5442(80) 


0.5411(103) 


Jo 3 


0.0529(10) 


0.0532(11) 



TABLE I: Comparison of observables measured during RHMC simulations with different specula- 
tive lower bounds for the spectrum of the quadratic Dirac operator: a) For m = 0.035, m = 0.3, 
P = 5.1370. b) For m = 0.03, m = 0.3, (3 = 5.1290. 

The good agreement between the 'data' produced with the different values for the spec- 
ulative lower bound indicates that our initial choice of 1 x 10~ 4 was adequate. One notes 
that, despite the high statistics - 300,000 trajectories for each run, the error bars are sizable. 
This is because the chosen f3 values are very close to /3 C , the transition j3, in each case which 
serves to maximize the fluctuations. 

Finally we compare the values of two fluctuation quantities, the Binder cumulant B^ipip) 
and the chiral susceptibility \^ a t Pc with results obtained from our HMD(R) simulations 
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over a range of dt values for m = 0.035 and an intermediate value of fij, \xi = 0.2, in figure 121 
The leading errors in HMD(R) are expected to be 0(dt 2 ). The fits linear in dt 2 in these 
figures appear to support this expectation, but indicate that we would need to include an 
0(dt 4 ) term to fit over the whole range of measurements. We note that the finite dt errors 
in both these graphs are quite large, so that agreement between HMD(R) and RHMC is 
non-trivial. 

IV. TESTS OF RHMC FOR xQCD 

Since we have been simulating the thermodynamics of xQCD with 2 massless quark 
flavours on 16 3 x 8 and 24 3 x 8 lattices, we use this larger lattice to run our tests. On this 
size lattice, extensive comparisons with the HMD(R) algorithm over a range of dt values is 
impractical, as is performing high statistics runs with various speculative lower bounds. 

Our first tests used a single-precision (32-bit floating-point) version of the code, then a 
modified version where the multimass conjugate gradient inversions which initiate and end 
each trajectory were replaced with double-precision (64-bit floating-point) versions, while 
leaving the inversion routines which perform the updates within each trajectory in single 
precision. We then performed reversibility tests of the same nature as those described in the 
previous section. For the updating, our convergence criterion was applied to the normalized 
conjugate gradient residual for the lowest lying pole in the partial fraction expansion of the 
rational approximation, while at the beginning and end of each trajectory it was applied to 
the unnormalized residual. (If Xq is an approximate solution to the linear equations Ax = b, 
the unnormalized residual is r u = \b — Axq\. The normalized residual is r n = r u /|6|.) We 
denote this upper limit on these residuals by r and note that, since our sources are extensive, 
this means that we are using a much more stringent convergence criterion at the beginning 
and end of each trajectory than in the updating, as required. 

Figure shows the changes in energy over the 'closed' trajectories t —>■ t + 1 —>■ t for 
the single precision code and for the same code with the initial and final inversions for each 
trajectory performed in double precision, r = 0.001 for these runs, since we have determined 
that no improvement is obtained by decreasing r below this value. (3 = 5.535, used for these 
simulations, is close to the critical point. The trajectory length is 1 time unit and dt = 0.025 
which gives an acceptance of ~ 75%. However, we turn off the accept/reject step in the 
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FIG. 2: a) The Binder cumulant B^tpip) as a function of dt 2 for the HMD(R) algorithm (crosses) 
and the RHMC algorithm (circle at dt = 0). b) The chiral susceptibility Xify as a function of dt 2 
for the HMD(R) algorithm (crosses) and the RHMC algorithm (circle at dt = 0). 
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RHMC algorithm for the purpose of these reversibility tests to speed up the passage of the 
system through phase space. The single precision code produces energy changes as large as 
~ 6 x 10~ 3 . More troubling is the fact that 5E is strongly biased to positive values. In fact the 
average SE is ~ 2.6 x 10 -3 which suggests that this code could produce small but significant 
systematic errors. The mixed single/double precision code shows more promise and we have 
extended our run to include 500 closed trajectories as shown in figure Et>. The maximum 
value of \SE\ is ~ 6.4 x 10~ 3 , which is larger than we would prefer. More problematic is 
that the average SE over the range is ~ 1.1 x 10 -3 , which suggests a systematic bias as 
is apparent in figure Eb- Unfortunately, we are limited by precision from improving this 
situation. We thus feel that to be safe we should increase the precision, especially if we wish 
to allow ourselves the option of increasing the lattice size in the future. 

We have tested reversibility of our double precision code as a function of the residual r. 
The results for 50 closed trajectories are shown in figure HJ We see that this measure of 
reversibility of the trajectories improves considerably as r is decreased from 0.01 to 0.001. 
Comparing this graph with that for the single precision tests, we see that precision is indeed 
the limiting factor. Further tests indicate that it is the convergence of the approximate 
inversions during the updating (rather than that of the more precise inversions at the start 
and end of each trajectory) that is the limiting factor. Careful monitoring of the trajectories 
for the r = 0.01 run, indicates that this is largely due to the fact that the number of 
iterations of the conjugate gradient inverter sometimes differed by 1 or 2 iterations between 
the corresponding steps of the forward and reversed trajectories, and this of course has the 
largest effect when one is far from convergence. This can be corrected by running for a fixed 
number of conjugate gradient iterations rather than at fixed r. Since at r = 0.01 it requires 
around 250 iterations on average for the inversion at each update, we performed a run of 1000 
trajectories with the number of inverter iterations at each update set to 300 (those at the 
ends of the trajectories were still set by a conservative choice of r). The results of running for 
1000 trajectories (500 closed) with the number of conjugate gradient iterations fixed at 300 
are shown in figure along with an identical number of trajectories with these inversions 
set by r = 0.002. The maximum value of \8E\ for the case where the number of conjugate 
gradient iterations is fixed at 300 is ~ 5 x 10~ 3 , no better than the single precision results. 
For r = 0.002 the maximum value of \5E\ is ~ 1.5 x 10~ 3 which we consider acceptable. The 
reason that fixing the number of conjugate gradient iterations at a relatively small value 
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XQCD 24 J x8 lattice, /3 = 5.535, 7=10 
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FIG. 3: a) Change in energy 5E for 50 closed trajectories for single precision code (single-single) and 
single precision updating with double precision initialization and measurement for each trajectory 
(single-double), b) 500 closed trajectories for single-double code. 
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XQCD 24 6 x8 lattice, £ = 5.535, y=10 
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FIG. 4: The change in energy 5E for closed trajectories for a range of choices of r for 64-bit 
precision code. 

did not give better results is presumably because this gives relatively poor convergence for 
some configurations. Studies of the HMC algorithm have indicated that if the convergence 
is too poor, the updating procedure develops instabilities which amplify the effects of finite 
precision errors 19]. We believe that such instabilities are most likely the reason why using 
a fixed number of conjugate gradient iterations during updating did not lead to smaller 
violations of reversibility. Figure shows the change of energies for the open trajectories 
from the same r = 0.002 run. 

We now need to justify our choices of spectral bounds on the quadratic Dirac operator. 
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XQCD 24 d x8 lattice, £ = 5.535, 7=10 
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FIG. 5: a) The change in energy 5E for closed trajectories for a fixed number of conjugate gradient 
iterations and for r = 0.002. b) The change in energy AE for open trajectories for r = 0.002. 

The chosen upper bound of 25 is obtained assuming that the magnitude of the average of 
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the 16 chiral auxiliary fields a + m which couple to the fermion bilinear ipif) at each site is 
bounded above by 1. Figure El shows the maximum over the lattice sites of the magnitude 
of a + m at the beginning and end of 3000 trajectories of our production running at m = 0, 
(3 = 5.535, 7 = 10, dt = 0.03125, trajectory length 1 and r = 0.002, where the acceptance is 
around 60-70%. As can be seen, 1 is a conservative upper bound for | a + Z7r|. 
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XQCD 24 J x8 lattice, £ = 5.535, y=10 
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FIG. 6: Maximum values of \a + in\ at the start and end of each trajectory of a production run. 



Since we are using a relatively large lattice (24 3 x 8) it is not practical to obtain adequate 
statistics close to the phase transition with more than one choice of speculative lower bounds 
as was done with our simulations of lattice QCD at finite isospin chemical potential, in order 
to check for consistent results. Instead, we performed a simulation of 1000 trajectories with 
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a speculative lower bound of 1 x 10 8 . At the beginning of each trajectory where we calculate 
the momenta conjugate to the pseudo-fermion fields, as 



we compared the values of p^ obtained with rational approximations which assume lower 
bounds of 1 x 10~ 6 , 1 x 1(T 8 and 1 x 1(T 10 . These used (25, 25), (30, 30), and (35, 35) rational 
approximations respectively, with maximum relative errors of 1.4 x 10 -11 , 3.0 x 10~ n and 
4.9 x 10~ n respectively. Similar comparisons were made of the £s at the end of each trajectory 
reconstituted as 



No comparisons with results from using different speculative lower bounds were made during 
the updates, since using different rational approximations corresponds to making different 
choices of the updating Hamiltonian, and such choices are unimportant unless they signifi- 
cantly affect the acceptance. 

In figure we show the magnitude of the difference between p^ values calculated using 
a lower bound of 1 x 10 -6 and those using a lower bound of 1 x 10~ 8 . We plot the maximum 
over sites and colours of this quantity for each trajectory as well as the average over sites 
and colours. Figure [7jp shows such differences for £. For 997 of these trajectories, the 
differences in p^, and £ values are consistent with zero within the accuracies of the rational 
approximations we use. For the other 3 trajectories, the differences are much larger. In 
addition, we note that the configurations which give relatively large errors for p^, and for 
£, are identical. Figure |H1 shows differences between the same quantities for lower bounds of 
10~ 8 and 10 -10 respectively for the same trajectories as in figure[7| In this case all differences 
are consistent with zero within the accuracy of the rational approximations we use. 

We interpret these results as indicating that for 997 of the configurations generated, all 
eigenvalues of the quadratic Dirac operator are > 10 -6 . For the remaining 3 configurations, 
the minimum eigenvalues lie between 10~ 8 and 10 -6 . It is reassuring to note that even for 
these configurations the differences are small, since a typical p^ or x has a magnitude 0(1). 
Hence it is probably reasonable to use a speculative lower bound of 10~ 6 and, unless this 
run was anomalous, choosing a speculative lower bound of 10~ 8 is completely safe. 

Finally we present some results from small lattice (8 3 x 4) simulations using the RHMC 
algorithm for xQCD and compare them with earlier published results using the HMD(R) 



(11) 



(12) 
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XQCD 24 3 x8 lattice, 10=5.535, 7=10 
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XQCD 24 3 x8 lattice, /S = 5.535, 7=10 
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FIG. 7: a) Magnitude of the difference between values using speculative lower bounds of 10 
and 10~ 8 . b) Magnitude of the difference between £ values using speculative lower bounds of 10 
and 10" 8 . 
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XQCD 24 d x8 lattice, 18=5.535, 7=10 
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and 10 -8 . b) Magnitude of the difference between £ values using speculative lower bounds of 10~ 10 
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algorithm Q|. For such small lattices precision is not an issue and these simulations were 
performed using 32-bit floating-point arithmetic. Table HTk gives the results for 7 = 10, 
m = for our RHMC simulations over a range of (3 values which bracket the phase transition. 
We ran 100,000 trajectories for each (3. Table ITTb gives the published results for the same 
parameters for the HMD(R) algorithm. Considering the relatively low statistics of the earlier 
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TABLE II: b) Observables measured in HMD(R) simulations on an 8 3 x 4 with 7 = 10, m = 0. 



HMD(R) results, the agreement with the RHMC simulations is excellent and the finite dt 
errors are small in this case. 



V. DISCUSSIONS AND CONCLUSIONS 

We have implemented the RHMC for 2 lattice actions where the lower bound on the 
spectrum of the quadratic Dirac operator is unknown. The first is the action for staggered 
lattice QCD with a finite chemical potential Hi for isospin, in the small fij regime. The 
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second is the xQCD action where the standard staggered quark action is augmented with 
a chiral 4-fermion interaction which allows simulations at zero quark mass. In the first case 
we simulate 3 fermion flavours, in the second 2. Both situations require taking a fractional 
power of the fermion determinant so that standard HMC methods fail, and the alternative 
HMD(R) algorithm gives finite dt errors. In each case the RHMC algorithm is implemented 
using a speculative lower bound to this spectrum, and its validity is tested by simulations. 
This rather naive approach appears successful in both cases. 

In each case we test the reversibility of our implementation. This would be exact for infi- 
nite precision arithmetic and exact inversions at the beginning and ending of the trajectories. 
We test this by reversing our trajectories and determining how close we come to the initial 
state. The measure of this is the difference between the classical energy of this final config- 
uration and that of the initial configuration. Since energy differences appear in the global 
Metropolis accept/reject step at the end of each trajectory, this energy difference should be 
much less than one. Because this energy is an extensive quantity, this requirement becomes 
more stringent as the lattice size is increased. For QCD at finite /xj, we are using 12 3 x 4 
lattices for which we find single-precision (32-bit) floating-point arithmetic to be adequate. 
For xQCD, where we are simulating on 24 3 x 8 lattices double-precision arithmetic appears 
desirable. Here we have also studied convergence requirements on the inversions needed to 
implement the rational approximations, both at the beginning and end of the trajectory and 
in the updating. Those at the beginning and end of each trajectory need to be performed 
with high precision. The inversions in the updating can be performed with less precision. If 
these updating inversions are too imprecise, the change in energy will be large and the final 
configuration will not be accepted. It appears, however, that the instabilities due to finite 
precision are worse when these inversions are too imprecise, and this is more important in 
determining the minimum acceptable precision than are acceptance criteria. 

For QCD at finite \xi we test our speculative lower spectral bounds by 2 methods. First 
we find that the RHMC results for Binder cumulants and chiral susceptibilities agree with 
extrapolations of our HMD(R) results to dt — 0, and that the position of the transitions 
calculated from both types of simulation appear consistent. Second, at the largest fij we 
use, we check that all observables and fluctuation quantities at the transition agree with 
those obtained with simulations using speculative lower bounds which are 10 or 100 times 
smaller than our initial choice. This indicates that a lower bound of 10~ 4 is adequate. 
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For xQCD, we find good agreement between measurements made using RHMC simula- 
tions and those made using HMD(R) simulations on 8 3 x 4 lattices. With our production 
lattice sizes of 24 3 x 8, it is impractical to make such detailed comparisons between the 2 
methods. In addition making runs with sufficient statistics to compare different speculative 
lower bounds is also too expensive. Instead, we have compared applying rational approxima- 
tions which assume different lower bounds to a set of 1000 configurations generated during 
one of our runs. We find agreement to within the accuracy of our rational approximations 
for speculative lower bounds of 10~ 8 and 10 -10 . An approximation using a speculative lower 
bound of 10~ 6 agrees with that using a speculative lower bound of 10~ 8 on most of these 
configurations. Even in those cases where the two estimates do not agree, the differences are 
small enough that the errors introduced in a production run would probably be negligible. 
We conclude therefore that a lower bound of 10~ 8 is almost certainly acceptable, and one 
of 10~ 6 is probably acceptable. We did apply this test to our QCD at finite pLj simulations 
comparing speculative lower bounds of 10~ 4 and 10~ 6 . Although there is no indication of 
any difference within the precision of our measurements, the limitation of 32-bit precision 
prevents us from using this test to make definitive statements. 

We conclude that the use of speculative lower bounds for the spectrum of the quadratic 
Dirac operator is adequate to allow use of the RHMC exact algorithm to simulate either of 
the theories we have considered. For theories where it is necessary to determine the lower 
bound dynamically, comparisons between the results of different speculative lower bounds 
could be made from a set of choices starting with that with the largest lower bound, during 
a run (as we have done), and choosing the first approximation which passes the test. This 
way the expensive calculations of the rational approximations using the Remez method, 
need only be performed once. The order of the rational approximation required to achieve 
the desired accuracy with a given lower bound does not increase very rapidly as that bound 
is increased. For example, in our applications where the relative error is always less than 
5 x 10 -11 , we achieve this for the 2-flavour case using a (20,20) rational approximation for 
a lower bound of 10~ 4 , a (25, 25) approximation for 10~ 6 , a (30, 30) approximation for 10~ 8 , 
and a (35, 35) approximation for 10~ 10 . Hence the 'library' of rational approximations one 
would need to maintain for such an implementation would be rather small. Such methods 
could be useful for domain-wall fermions and overlap fermions which have the problem that 
their lower spectral bounds are also unknown, and might need to be determined dynamically 
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